Data on tour

Last week I was on Gran Canaria for a vacation We took hiking trips and explored the island by car and on foot. So, what better way to keep up the holiday spirit a while longer than to visualize all the places we went!?

I am combining the location data collected by

our car GPS,
Google location data from my phone and
the hiking tracks we followed.

The map

I am using ggplot and ggmap. I centered the map to the midpoint of Gran Canaria as given by http://www.travelmath.com/island/Gran+Canaria, which is 27° 57’ 54" N / 15° 35’ 59" W. Converting to decimal with http://www.rapidtables.com/convert/number/degrees-minutes-seconds-to-degrees.htm that is 27.965 & 15.59972.

library(ggplot2)
library(ggmap)

map_theme <- list(theme(legend.position = "top",
                        panel.grid.minor = element_blank(),
                        panel.grid.major = element_blank(),
                        panel.background = element_blank(),
                        plot.background = element_rect(fill = "white"),
                        panel.border = element_blank(),
                        axis.line = element_blank(),
                        axis.text.x = element_blank(),
                        axis.text.y = element_blank(),
                        axis.ticks = element_blank(),
                        axis.title.x = element_blank(),
                        axis.title.y = element_blank(),
                        plot.title = element_text(size = 18)))

map <- get_map(c(lon = -15.59972, lat = 27.965), zoom = 10, maptype = "terrain-background", source = "google")
#map <- get_map(c(lon = -15.59972, lat = 27.965), zoom = 10)

GPS tracks

The GPS tracks in .gpx format were downloaded from the device (Garmin) the same way as I had already done with the tracks from our US/Canada roadtrip last year. Garmin’s .gpx tracks were not in a standard format that could be read with readGPX(), so I had to use htmlTreeParse() from the XML library.

library(XML)

myfiles <- list.files(path = "Takeout_2017_March", full.names = TRUE)
myfiles <- myfiles[grep(".gpx", myfiles)]

for (i in 1:length(myfiles)) {
  tryCatch({
    pfile <- htmlTreeParse(myfiles[i], useInternalNodes = T)

  }, error = function(e){cat("ERROR\n")})

  if (exists("pfile")) {
    # Get all elevations, times and coordinates via the respective xpath
    elevations <- as.numeric(as.character(xpathSApply(pfile, path = "//trkpt/ele", xmlValue)))
    times <- xpathSApply(pfile, path = "//trkpt/time", xmlValue)
    coords <- xpathSApply(pfile, path = "//trkpt", xmlAttrs)
    # Extract latitude and longitude from the coordinates
    lats <- as.numeric(as.character(coords["lat",]))
    lons <- as.numeric(as.character(coords["lon",]))
    # Put everything in a dataframe and get rid of old variables
    geodf <- data.frame(lat = lats, lon = lons, ele = elevations, time = times)

    if (i == 1) {
      geodata <- geodf
      } else {
      geodata <- rbind(geodata, geodf)
    }
  }
  rm(pfile)
}

geodata_all <- geodata

# Transforming the time column
geodata_all$time <- as.character(strptime(geodata_all$time, format = "%Y-%m-%dT%H:%M:%SZ"))

# ordering by date
library(dplyr)
geodata_all <- arrange(geodata_all, time)

geodata_gc <- filter(geodata_all, lat < 28.5 & lat > 27.5) %>%
  filter(lon > -16 & lon < -15)

GPX hiking tracks

The hiking tracks we followed came mostly from the Rother Wanderführer, 7th edition from 2016. They were in standard .gpx format and could be read with readGPX().

library(plotKML)

myfiles <- list.files(path = "hiking", full.names = TRUE)
myfiles <- myfiles[grep(".gpx", myfiles)]

for (i in 1:length(myfiles)) {
  t <- readGPX(myfiles[i])
  
  geodf <- data.frame(lon = t$tracks[[1]][[1]]$lon,
                      lat = t$tracks[[1]][[1]]$lat,
                      ele = t$tracks[[1]][[1]]$ele,
                      tour = names(t$tracks[[1]]))

    if (i == 1) {
      geodata <- geodf
      } else {
      geodata <- rbind(geodata, geodf)
    }
  rm(pfile)
}

geodata_tracks <- geodata

Google location history

The Google location data I loaded the same way as in my post about plotting your Google location history.

library(jsonlite)
system.time(x <- fromJSON("Standortverlauf.json"))

# extracting the locations dataframe
loc = x$locations

# converting time column from posix milliseconds into a readable time scale
loc$time = as.POSIXct(as.numeric(x$locations$timestampMs)/1000, origin = "1970-01-01")

# converting longitude and latitude from E7 to GPS coordinates
loc$lat = loc$latitudeE7 / 1e7
loc$lon = loc$longitudeE7 / 1e7

# keep only data from Gran Canaria
loc_gc <- filter(loc, lat < 28.5 & lat > 27.5) %>%
  filter(lon > -16 & lon < -15)

All three datasets had information about the latitude/longitude coordinate pairs and of the elevation of each observation, so I can combine these three attributes from the three location data sources.

geodata_tracks$ele <- as.numeric(as.character(geodata_tracks$ele))
data_combined <- data.frame(lat = c(geodata_gc$lat, loc_gc$lat, geodata_tracks$lat),
                            lon = c(geodata_gc$lon, loc_gc$lon, geodata_tracks$lon),
                            ele = c(geodata_gc$ele, loc_gc$altitude, geodata_tracks$ele),
                            track = c(rep("GPS", nrow(geodata_gc)), rep("Google", nrow(loc_gc)), rep("Hiking", nrow(geodata_tracks))))
data_combined <- data_combined[!duplicated(data_combined), ]

ggmap(map) + 
  geom_point(data = data_combined, aes(x = lon, y = lat, color = track), alpha = 0.3) +
  scale_color_brewer(palette = "Set1") +
  map_theme

ggmap(map) + 
  geom_point(data = na.omit(data_combined), aes(x = lon, y = lat, color = ele)) +
  scale_color_gradient2(low = "lightblue", mid = "blue", high = "red", name = "Elevation") +
  map_theme

library(dplyr)
geodata_gc$time <- as.POSIXct(geodata_gc$time, format = "%Y-%m-%d %H:%M:%S", tz = "GMT")

data_combined_2 <- geodata_gc

# removing duplicate entries
data_combined_2 <- data_combined_2[!duplicated(data_combined_2), ]

# order by time
data_combined_2 <- data_combined_2[order(data_combined_2$time, decreasing = FALSE), ]

# Shifting vectors for latitude and longitude to include end position
shift.vec <- function(vec, shift) {
  if (length(vec) <= abs(shift)) {
    rep(NA ,length(vec))
  } else {
    if (shift >= 0) {
      c(rep(NA, shift), vec[1:(length(vec) - shift)]) }
    else {
      c(vec[(abs(shift) + 1):length(vec)], rep(NA, abs(shift)))
    }
  }
}

data_combined_2$lat.p1 <- shift.vec(data_combined_2$lat, -1)
data_combined_2$lon.p1 <- shift.vec(data_combined_2$lon, -1)

# Calculating distances between points (in metres) with the function pointDistance from the 'raster' package.

library(raster)
data_combined_2$dist.to.prev <- apply(data_combined_2, 1, FUN = function(row) {
  pointDistance(c(as.numeric(as.character(row["lat.p1"])),
                  as.numeric(as.character(row["lon.p1"]))),
                c(as.numeric(as.character(row["lat"])), as.numeric(as.character(row["lon"]))),
                lonlat = T) # Parameter 'lonlat' has to be TRUE!
})

round(sum(as.numeric(as.character(data_combined_2$dist.to.prev)), na.rm = TRUE) * 0.001, digits = 2)

## [1] 1009.6

library(maptools)

## Loading required package: sp

## Checking rgeos availability: TRUE

library(raster)
srtm <- getData("SRTM", lon=-15.59972, lat=27.965)

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

raster::persp(srtm,
              xlab = "Longitude", ylab = "Latitude", zlab = "Elevation",
              phi = 45, theta = 0,
              expand = 0.1,
              col = "forestgreen", border = NULL,
              box = TRUE)

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

e <- extent(min(data_combined$lon) - 0.2, # xmin
            max(data_combined$lon) + 0.1, # xmax
            min(data_combined$lat) - 0.1, # ymin
            max(data_combined$lat) + 0.1) # ymax

srtm_c <- crop(srtm, e)

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

raster::persp(srtm_c,
              xlab = "Longitude", ylab = "Latitude", zlab = "Elevation",
              phi = 45, theta = 45,
              expand = 0.5,
              col = "forestgreen", border = NULL,
              box = TRUE)

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

raster::persp(srtm_c,
              xlab = "Longitude", ylab = "Latitude", zlab = "Elevation",
              phi = 45, theta = 225,
              expand = 0.5,
              col = "forestgreen", border = NULL,
              box = TRUE)

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

x <- terrain(srtm_c, opt = c("slope", "aspect"), unit = "degrees")

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

plot(x)

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

slope <- terrain(srtm_c, opt = "slope")

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

aspect <- terrain(srtm_c, opt = "aspect")

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

hill <- hillShade(slope, aspect, 40, 270)

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

plot(hill, col = grey(0:100/100))

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

library(rgdal)
library(rasterVis)
library(rgl)
library(htmlwidgets)

options(rgl.printRglwidget = TRUE)

open3d()

## glX 
##   1

plot3D(srtm)

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

open3d()

## glX 
##   2

plot3D(srtm_c)

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

rglscene <- scene3d()
# render the saved rgl scene as widget in the markdown page
rglwidget(rglscene)

rgl.close()

Data on tour

Dr. Shirin Glander

2017-03-28

The map

GPS tracks

GPX hiking tracks

Google location history